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Abstract 

This  report  demonstrates  the  snow  water  equivalent  (SWE)  modeling 
capabilities  present  within  the  Gridded  Surface  Sub-surface  Hydrologic 
Analysis  (GSSHA)  model  (Downer  and  Ogden  2004,  2006)  by  comparing 
simulated  and  observed  snow  accumulation  data.  GSSHA  is  widely  used  in 
both  military  and  civil  works  projects  to  simulate  runoff,  streamflow,  and 
sediment  and  constituent  transport.  GSSHA  originally  employed  an  energy 
balance  method  for  calculating  snow  accumulation  and  melt,  but  now  has 
an  optional  temperature-index  method  and  hybrid  energy  balance  method 
(explained  within  this  report)  as  well.  The  full  snowpack  modeling 
capabilities  of  the  Snowmelt  Numerical-Analytical  Package  (SNAP)  model 
(Albert  and  Krajeski  1998)  have  also  been  included  to  simulate  snow  pack 
depth  and  density  regardless  of  which  of  the  three  snow  accumulation 
models  are  employed.  With  proficient  SWE  modeling  capabilities,  the 
GSSHA  model  can  be  utilized  to  solve  complex  watershed-related  issues  in 
regions  where  snow  accumulation  and  melt  are  often  the  most  critical 
source  of  water. 
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Preface 

This  report  describes  the  procedures  and  results  from  a  study  to  improve 
the  accuracy  of  snow  accumulation  and  melting  process  simulation  within 
the  Gridded  Surface  Subsurface  Hydrologic  Analysis  (GSSHA)  model.  The 
study  was  conducted  from  October  2011  to  September  2012. 

This  study  was  conducted  for  the  Flood  and  Coastal  Storm  Damage 
Reduction  Research  and  Development  Program  under  “Enhancing  the 
Snowmelt  Routines  within  GSSHA.”  The  Program  Manager  was 
Dr.  Cary  A.  Talbot. 

The  work  was  performed  by  Michael  L.  Folium  and  Dr.  Charles  W.  Downer, 
Hydrologic  Systems  Branch,  (HSB),  Flood  and  Storm  Protection  Division 
(FSPB),  US  Army  Engineer  Research  and  Development  Center  (ERDC)  - 
Coastal  and  Hydraulics  Laboratory  (CHL). 

Dr.  Aaron  R.  Byrd  was  Chief,  HSB;  Dr.  Ty  V.  Wamsley  was  Chief,  FSPD;  and 
William  R.  Curtis  was  Technical  Director,  Flood  and  Coastal  Storm  Damage 
Reduction  Research  and  Development  Program.  Jose  E.  Sanchez  was  Acting 
Director,  CHL. 

At  the  time  of  publication,  COL  Kevin  J.  Wilson  was  Commander  of  ERDC, 
and  Dr.  Jeffery  P.  Holland  was  Director. 
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1  Introduction 

The  US  Army  Engineer  Research  and  Development  Center  (ERDC)  has  a 
long  history  in  developing  and  applying  gridded  physics-based  models  to 
help  solve  complex  watershed-related  issues.  ERDC’s  Coastal  and 
Hydraulics  Laboratory  (CHL)  developed,  and  continues  to  further  enhance, 
the  capabilities  of  the  continuous  simulation,  physics-based,  spatially 
distributed,  hydrologic  model  GSSHA.  A  simplistic  energy  balance  method 
was  originally  employed  by  GSSHA  for  determining  snow  accumulation  and 
melt.  Applications  of  this  method  demonstrated  both  volume  and  timing 
errors  for  snow  melt.  This  report  outlines  recent  improvements  made  to  the 
snow  accumulation  and  melt  routines  within  GSSHA  to  make  it  more 
applicable  in  snow-dominated  regions. 

Three  methods  are  described  in  this  report;  original  energy  balance, 
hybrid  energy  balance,  and  temperature  index.  These  are  available  within 
the  GSSHA  model  to  simulate  snow  accumulation  and  melt.  By  way  of 
examples,  how  well  these  methods  perform  in  terms  of  comparisons  with 
observed  data  at  several  test  sites  in  the  western  United  States  is  shown. 
This  report  also  shows  how  accounting  for  the  temperature  of  the 
snowpack,  known  as  heat  deficit  or  cold  content,  can  significantly  improve 
the  accuracy  of  the  simulations. 

1.1  GSSHA  model  description 

GSSHA  is  a  spatially-distributed,  physics-based,  continuous  simulation 
hydrologic  model.  Computations  in  GSSHA  are  performed  on  a  uniform 
raster  grid,  which  is  used  to  represent  the  modeled  watershed  system 
(Figure  l).  Numerous  hydrologic  components  -  such  as  rainfall  distribution, 
evapotranspiration  (ET),  infiltration,  and  snow  melt/accumulation  -  are 
calculated  within  each  grid  cell  of  the  model  at  a  user-defined  time  step. 
Although  these  hydrologic  components  are  calculated  within  each 
individual  cell,  they  are  connected  through  two-dimensional  (2D)  overland 
flow,  2D  groundwater  flow,  and  a  one-dimensional  stream  network,  which 
results  in  the  overall  watershed  response.  GSSHA  also  has  the  capacity  to 
simulate  interaction  between  surface  water  and  groundwater  through 
infiltration  and  connections  between  the  stream  network  and  the 
groundwater. 
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Figure  1.  Structured  grid  format  of  the  GSSHA  Trout  Creek  river  basin,  California  model 

(elevation  in  meters). 


Hydrometeorology  (HMET)  data,  as  well  as  soil  and  land  use  charac¬ 
teristics,  are  unique  to  each  individual  grid  cell,  making  the  model  useful  in 
situations  where  uniform  assumptions  (i.e.,  the  entire  basin  has  the  same 
land  use  characteristics,  soil  type  characteristics,  precipitation,  etc.)  are  not 
realistic.  Altering  the  snow  accumulation  and  melt  algorithms  will  change 
the  way  they  are  calculated  within  each  grid  cell,  affecting  them  for  the 
entire  watershed. 

In  the  original  GSSHA  snow  accumulation  and  melt  model,  precipitation 
was  assumed  to  be  frozen  when  the  air  temperature  was  below  o°C. 
Melting  of  the  snow  was  assumed  to  occur  when  sufficient  energy  was 
transferred  to  the  snowpack  from  a  variety  of  sources:  precipitation, 
sensible  heat  transfer  due  to  turbulence,  evaporation  and  sublimation, 
atmospheric  longwave  radiation,  and  longwave  radiation  emission  by  the 
soil.  Because  the  energy  balance  formulation  did  not  account  for  the  sub¬ 
freezing  temperatures  within  the  snowpack,  it  likely  overestimated  melt, 
therefore,  underestimating  the  snowpack.  Although  the  surface  tempera¬ 
ture  might  rise  above  o  °C,  melting  within  the  pack  will  not  occur  until  the 
temperature  of  the  snowpack  rises  above  o  °C.  The  heat  that  must  be 
applied  to  bring  the  snowpack  temperature  to  o°C  is  referred  to  as  the 
heat  deficit. 
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1.2  Heat  deficit 

The  heat  deficit,  also  referred  to  as  cold  content,  is  defined  as  the  amount 
of  heat  that  must  be  added  to  return  the  snowpack  from  below  o  °C  to  an 
isothermal  state  (o  °C)  (Melloh  1999;  Anderson  2006).  The  equations  to 
calculate  the  heat  deficit  within  this  report  are  based  on  the  National 
Weather  Service  River  Forecasting  System  (NWSRFS)  SNOW-17  model 
(Anderson  1973,  2006).  SNOW-17  employs  two  factors  that  can  increase  or 
decrease  the  heat  deficit:  (1)  energy  flux  between  the  snowpack  and  its 
surroundings,  such  as  a  temperature  gradient  between  the  atmospheric 
temperature  and  the  antecedent  temperature  of  the  snowpack,  and  (2) 
mass  change,  such  as  precipitation.  It  is  important  to  include  the  heat 
deficit  because  melt  and  rain  continue  to  refreeze  within  the  snow  cover 
until  the  heat  deficit  reaches  zero  (Melloh  1999).  The  heat  deficit  described 
below  is  in  units  of  millimeters  of  SWE,  making  it  easily  incorporated  into 
most  snow  melt  routines  by  simply  reducing  the  amount  of  melt  (in  mm  of 
SWE)  calculated  in  the  melt  routine  by  the  heat  deficit. 

1.2.1  Energy  flux  heat  deficit  due  to  temperature  gradients 

Temperature  gradients  between  the  atmosphere  and  the  snowpack  drive 
diffusion  of  heat  within  the  snowpack,  either  increasing  or  decreasing  the 
heat  deficit.  When  the  atmospheric  temperature  drops  below  the  antecedent 
temperature  of  the  snowpack,  the  snowpack  temperature  drops  as  well, 
increasing  the  temperature-caused  heat  deficit.  When  the  atmospheric 
temperature  rises  above  the  antecedent  temperature  of  the  snowpack,  the 
snow  pack  temperature  gradually  increases  until  it  reaches  o  °C,  its 
maximum  value.  This  increase  in  the  snowpack  temperature  decreases  the 
temperature-caused  heat  deficit  until  it  reaches  zero. 

Ignoring  the  temperature  deficit  of  the  snow  causes  the  simulated  snow  to 
melt  too  quickly  when  the  air  temperature  rises  from  below  o  °C  to  above 
o  °C.  The  SNOW-17  model  incorporates  a  method  for  which  temperature 
indices,  essentially  a  term  to  consider  the  snowpack  temperature,  are 
calculated  based  on  Equation  1,  and  then  are  used  in  Equation  2  to  calculate 
the  change  in  snow  cover  heat  deficit  due  to  temperature.  Equations  1,  2, 
and  3  are  used  to  account  for  the  temperature  deficit  within  the  snowpack. 
The  weighting  multiplier  (27PM)  and  the  proportionality  factor  ( NMf )  are 
calibration  parameters,  but  results  presented  later  in  this  report  will  show 
they  are  relatively  insensitive,  compared  to  the  parameters  and  algorithms 
pertaining  to  the  melting  processes. 
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ATI,  =ATI,_1  +TIPMdl  (Ta - ATI 

(1) 

AD,  =NMr{ATI,_1  -Ta) 

(2) 

TIPMj,  =  1.0  -  (1.0  -  TIPMj"" 

(3) 

where 

ATIi  =  temperature  indices  at  i  time  step  (°C);  at  initial  time  step 
ATIo  =  o.o 

Ta  =  air  temperature  (°C) 
dt  =  time  step  (hours) 

TIPMdt  =  weighting  multiplier  (o.i  to  l.o)  for  previous  time  period 
TIPM  =  weighting  multiplier  (o.i  to  l.o)  for  previous  6-hr  period 
A Dt  =  change  in  heat  deficit  due  to  a  temperature  gradient,  expressed 
in  water  equivalent  (mm  SWE) 

NMf  =  proportionality  factor,  also  known  as  the  negative  melt  factor 
(mm/°C  per  6  hr) 

1.2.2  Mass  change  heat  deficit  due  to  precipitation 

Like  temperature  discussed  above,  precipitation  can  also  affect  the  heat 
deficit.  Precipitation  can  increase  or  decrease  the  mass  of  the  snowpack 
and  depending  on  the  temperature  of  the  precipitation,  can  increase  or 
decrease  the  heat  deficit.  Equation  4  is  used  to  capture  the  precipitation 
portion  of  the  heat  deficit  within  the  snow  cover  while  precipitation  is 
occurring.  Continuous  accounting  of  the  heat  deficit  due  to  both 
temperature  and  precipitation  is  computed  with  Equation  5. 


ADp=(p.*r,)/(t/ /<=<.) 

(4) 

A  =  A  1 +AA  +^A> 

(5) 

where 

A Dp  =  change  in  heat  deficit  due  to  precipitation,  expressed  in  water 
equivalent  (mm  SWE) 

Pn  =  water  equivalent  of  adjusted  precipitation  (mm)  (see 
Section  2.1) 
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TP 

Lf 

?ice 

Di 


temperature  of  snowfall;  assumed  equal  to  air  temperature  or 
o  °C,  whichever  is  less  (°C) 
latent  heat  of  fusion  (80  caPgnr1) 
specific  heat  of  ice  (0.5  caPgnr1*  °C1) 

heat  deficit  of  snowpack  during  i  time  step,  expressed  in  water 
equivalent  (mm  SWE);  at  initial  time  step  D0- 0.0 
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2  Description  of  the  Snow  Accumulation 
and  Melt  Methods 

All  three  methods  used  within  GSSHA  to  simulate  snow  melt  (energy 
balance,  hybrid  energy  balance,  and  temperature-index)  assume  the 
snowpack  consists  of  a  single  layer.  Certain  capabilities  such  as  time 
variations  of  liquid  water  content  (Bloschl  and  Kirnbauer  1991),  interflow 
within  the  snowpack  layers  due  to  ice  sheets,  and  avalanche  modeling 
(Colbeck  1991)  do  exist  when  applying  multi-layer  snowpack  models;  but 
the  required  data  to  apply  such  models  at  the  watershed  scale  is  unrealistic 
in  most  basins.  Multi-layer  snow  models  are  typically  deployed  at  the  site- 
scale  where  spatially  close  data  are  available,  making  them  impractical  in 
many  sparse-data  cases. 

In  nature,  snow  accumulation  and  melt  are  part  of  a  distributed  process. 
Due  to  the  spatially  distributed  structure  of  the  GSSHA  model,  the  snow  is 
modeled  as  a  distributed  process  that  gives  it  more  potential  for 
addressing  a  variety  of  real-world  problems  rather  than  a  semi-distributed 
or  lumped-parameter  hydrologic  model  (Kirnbauer  et  al.  1994). 

2.1  Snow  accumulation 

Independent  of  the  three  melt  methods,  GSSHA  has  a  means  to  account 
for  inaccuracies  in  the  gauge  systems  as  well  as  the  temperature  at  which 
snow  begins  to  form.  Typically,  when  the  air  temperature  is  at  or  below 
o°C,  precipitation  falls  in  the  form  of  snow.  However,  due  to  temperature 
differences  at  snow  formation  height  (i.e.  temperature  in  the  upper 
atmosphere  where  snow  forms  may  be  lower  than  the  ground  tempera¬ 
ture),  snow  can  form  when  the  air  surface  temperatures  are  greater  than 
o°C.  For  this  reason,  the  temperature  at  which  precipitation  begins  to  fall 
as  snow  ( MBASE ,  o°C)  is  considered  a  global  calibration  parameter. 
During  precipitation,  when  the  air  temperature  is  below  MBASE,  land 
surface  accumulation  is  assumed  to  be  snow  or  ice.  If  snow  is  already 
present  in  a  cell,  the  accumulation  is  added  to  the  existing  snow. 

Due  to  measure  catch  deficiencies,  precipitation  gauges  often  underestimate 
the  amount  fallen  in  the  form  of  snow  (Anderson  2006).  As  a  result,  the 
snow  adjustment  factor  (SCF)  is  applied  to  the  gauge  measurement  of 
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precipitation  (Px)  to  create  the  adjusted  precipitation  of  the  newly  fallen 
snow  (Pn)  (see  Equation  6).  SCF is  considered  a  calibration  parameter, 
while  the  fraction  of  precipitation  in  the  form  of  snow  (fs )  is  considered 
constant  at  1.0  when  temperatures  are  at  or  below  MBASE,  and  o.o  when 
above  MBASE.  This  in  effect  causes  the  model  to  treat  partially  frozen 
precipitation  such  as  frozen  rain,  commonly  referred  to  as  “slush”  or  “sleet,” 
as  snow  or  rain,  depending  on  the  temperature  and  MBASE. 

By  not  deciphering  when  precipitation  is  in  the  form  of  frozen  rain,  the 
snowpack,  under  specialized  circumstances,  might  be  slightly  overestimated 
in  the  beginning  stages.  For  example,  if  frozen  rain  occurs  when  there  is  not 
a  snowpack,  the  model  could  simulate  the  frozen  rain  as  snow  and  form  a 
snowpack,  when  in  reality  none  are  formed.  If  a  snowpack  already  exists 
during  a  frozen  rain  event,  there  will  be  no  effect  on  the  snow  simulation 
because  the  temperature  and  volume  of  the  snowpack  and  precipitation  are 
the  only  factors  that  contribute  to  the  increase  or  decrease  of  the  snowpack; 
whether  the  precipitation  is  in  the  form  of  snow,  rain,  or  frozen  rain  makes 
no  difference  (Section  2.2).  Any  errors  associated  with  the  assumption  that 
precipitation  is  in  the  form  of  snow  or  rain  will  be  infrequent  and  small  due 
to  the  rarity  of  occurrences  (frozen  rain  simulated  as  snow  when  no 
previous  snowpack  exists)  and  relatively  small  volumes  of  precipitation 
during  these  occurrences. 


P„=P'f'SCF 


(6) 


where 

Px  =  gauge  measurement  of  water  equivalent  of  precipitation  (mm) 
fs  =  fraction  of  precipitation  in  the  form  of  snow 
SCF  =  snow  adjustment  factor  (unitless) 

2.2  Snow  melt 

2.2.1  Energy  balance  (EB)  method 

The  energy  balance  (EB)  method  of  estimating  snow  melt  uses  a  simple 
conceptualization  of  the  snowpack  (Downer  and  Ogden  2006).  This 
method  calculates  the  total  energy,  (cal  cnr2  hr1),  available  hourly  to  melt 
snow  based  on  energy  fluxes  due  to: 


1.  Evaporation  and  sublimation  (Qe,  Equations  7  and  8); 
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2.  Sensible  heat  transfer  due  to  turbulence  (Qh,  Equation  9); 

3.  Atmospheric  longwave  radiation  (Qa,  Equation  10); 

4.  Longwave  emission  by  soil  (constant,  Qbs  =  27.0);  and 

5.  Precipitation  (Qp,  Equation  11) 

*all  energy  fluxes,  (Q),  are  in  units  of:  cal  cm  2  hr1 


f(v)  =  0.0002  *  u 

(7) 

vh 

Qe=f(V)*"l*es  6.11 

e  v  ’  100 

(8) 

Qh=0.366*Ta*Pa*f(V) 

(9) 

,  =  49.56xlO“10  *(Ta  +273)4 

(10) 

p  *T 

Q  -  n  a 

100 

(11) 

where 

f[V)  =  wind  function  (km  hr1)  (Anderson  1978) 

U  =  wind  speed  (m  S’1) 
rh  =  relative  humidity  (%) 

6sat  =  saturated  vapor  pressure  (mb) 

Pa  =  atmospheric  pressure  (mb) 

Equation  12  calculates  the  total  energy  to  melt  snow  using  these  energy 
fluxes.  The  latent  heat  of  melting,  (~8o  cal  cm-3),  is  then  used  to  calculate 
the  melted  water  (Equation  13). 


'melt  ~Qa  ~Qbs  +Qe+Qh+Qp 

(12) 

MEB-  Qmelt/100^  *dt 

80.0 

(13) 

where 


Qmeit  =  total  energy  available  to  melt  snow  (cal  cm-2  hr1) 
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Meb  =  melt  calculated  using  energy  balance  snow  melt  routine 
(mm  SWE) 

The  EB  method  of  estimating  snowfall  accumulation  and  melting  does  not 
include  any  representation  of  the  snowpack’s  temperature.  This  method 
also  assumes  a  constant  soil  temperature,  resulting  in  a  constant  value  for 
longwave  emission  by  the  soil  (Qbs= 27.0  cal  cm-1  hr1).  Future  work  efforts 
might  look  at  calculating  this  value,  but  a  soil  temperature  model  is 
required.  More  information  on  the  methods  and  equations  (Equations  7 
through  13)  used  to  calculate  the  energy  fluxes  in  this  method  can  be  found 
on  the  GSSHA  Wiki  site  (http://www.gsshawiki.com). 


2.2.2  Temperature-index  (Tl)  method 

The  temperature-index  (TI)  method  of  estimating  snow  melt  is  based  on 
the  NWSRFS  SNOW-17  model.  This  method  takes  into  consideration  the 
time  of  year,  melt  due  to  temperature,  melt  due  to  precipitation,  and  heat 
deficits  within  the  snowpack.  Equations  14  and  15  are  used  to  calculate  the 
amount  of  melt  during  a  time  step.  Equation  14  is  used  in  precipitation- 
dominated  time  spans  when  the  average  precipitation  over  the  previous  6 
hours  has  exceeded  0.25  mm  hr1  and  precipitation  is  occurring  during  the 
current  time  step.  Equation  15  is  used  when  the  melt  is  considered 
temperature-dominated.  The  precipitation  temperature,  Tr,  is  assumed  to 
be  o  °C  or  the  temperature  of  the  air,  whichever  is  greater.  Both  melt 
routines  work  only  when  the  air  surface  temperature  is  greater  than  o°C 
and  only  one  melt  routine  is  run  per  time  step. 


This  method  also  keeps  track  of  melt  water  being  stored  and  water  being 
released  from  the  snowpack.  Although  Equations  1  through  5  help 
simulate  the  ripeness  of  the  snowpack  by  accounting  for  the  heat  deficit, 
Equations  18  through  21  are  included  in  the  TI  method  to  help  determine 
how  much  liquid  water  is  being  stored  within  the  snowpack.  The  amount 
of  melt  that  leaves  the  snowpack  (Mr/)  is  calculated  in  Equation  20  and  is 
representative  of  the  overall  melt  from  the  pack  after  the  heat  deficit  and 
water  storage  capacity  are  accounted  for. 


M„  =o*df 


(Ta+  273)4— 2734 


+  0.0125  *PX* 


A 


T  + 

use  r  1 


8.5  *  fua*  ( dt  /  6)  *  [(0.9  *  esat  -  6.1l)  +  0.00057  *Pa*Ta 


(14) 
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Mnr  =  Mf*  (Ta  -  MBASE )  +  0.0125  *  Px*  fruse  *  Tr 
Mf  =  (dt  /  6) *  [S„  *  A,  *  ( MFmax  -  MF„m )  +  MFm[n  } 


■  0.5  *  sin 


N*  2*n 


366 


+  0.5 


Wqx  —  PLWHC  *  Wi 

Qw  =  Mr  +  Mnr  +  Px*  fvuse 
Mti  =Qw+  wq,n  -  Wqx  -D2-  (PLWHC *  D2) 
Wq>n+1=Wqtn-Qw-D2 


(15) 

(16) 

(17) 

(18) 

(19) 

(20) 
(21) 


where 

Mr  = 
o  = 

fruse  = 

fra  = 
Mnr  = 

Mf  - 
MFmin  = 
MFmax  = 
MBASE  = 
Av  = 
Sv  = 
N  = 
Mti 

Wqx  = 
Wq  = 
Wi  = 
Qw  = 
PLWHC  = 


melt  during  precipitation-dominated  time  steps  (mm) 
Stefan-Boltzman  constant  (6.12*io-10  mm/°K/hr) 
fraction  of  precipitation  in  the  form  of  rain 
average  wind  function  (mm/mb  per  6  hr) 
melt  during  temperature-dominated  time  steps  (mm) 
melt  factor  (mm/°C/dt) 

minimum  melt  factor  on  -  Dec  21  (mm  SWE/°C  per  6  hr) 
maximum  melt  factor  on  -  June  21  (mm  SWE/°C  per  6  hr) 
temperature  at  which  precipitation  begins  to  fall  as  snow  (°C) 
seasonal  variation  adjustment 
seasonal  sine  curve  melt  variation 
day  number  since  March  21 

overall  melt  calculated  using  the  temperature-index  snow  melt 

routine  (mm  SWE) 

liquid  water  capacity  (mm) 

liquid  water  held  by  the  snow  (mm) 

water  equivalent  of  the  ice  portion  of  the  snow  cover  (mm) 
liquid  water  available  at  the  snow  surface  (mm) 
percent  liquid  water  holding  capacity  (decimal  fraction) 


The  ice  portion  of  the  snowpack  (Wi)  is  calculated  based  on  the  SWE  and 
PLWHC  of  the  snowpack.  The  calibration  parameters  for  the  TI  method 
include:  SCF,  MFmax,  MFmin,  fra,  MBASE,  fruse,  TIPM,  NMf,  and  PLWHC. 
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For  more  information  on  the  SNOW-iy  model  and  the  equations  presented 
(14  through  21),  the  reader  is  guided  to  Anderson  (1973  and  2006)  and 
Melloh  (1999).  The  SNOW-1 7  model  was  intended  for  a  6-hr  time  step,  but 
is  compatible  with  an  hourly  time  step.  Although  GSSHA  employs  a  global 
variable  time  step,  which  can  be  sub-minute,  the  TI  method  is  run  hourly. 
The  melt  generated  during  the  hour  time  step  is  distributed  to  the  other 
GSSHA  model  processes  at  the  global  variable  time  step. 

2.2.3  Hybrid  energy  balance  (HY)  method 

A  simple  adjustment  to  the  EB  method  was  made  to  create  the  hybrid 
energy  balance  (HY)  method  for  snow  melt.  By  accounting  for  the  heat 
deficit  within  the  melting  algorithms  of  the  EB  method,  the  melting 
processes  and  the  snow  pack  dynamics  (in  the  form  of  heat  deficit)  are 
accounted  for,  making  the  model  more  accurate.  Equation  22  shows  the 
simple  adjustment  that  creates  the  HY  method.  Equations  1  through  5  show 
how  the  heat  deficit  is  calculated.  The  HY  method  requires  few,  if  any, 
calibration  parameters  because  the  EB  method  does  not  have  any  and  the 
calibration  parameters  associated  with  the  heat  deficit  equations  are 
relatively  insensitive.  Although  GSSHA  uses  a  global  variable  time  step, 
which  can  be  sub-minute,  the  heat  deficit  is  calculated  at  an  hourly  time 
step.  Therefore,  the  HY  method  is  also  run  at  an  hourly  time  step.  The  melt 
generated  during  the  hourly  time  step  is  distributed  to  the  other  GSSHA 
model  processes  at  the  global  variable  time  step. 

Mhy=Meb-D2  (22) 


where 


Mhy  =  melt  calculated  using  hybrid  snow  melt  routine  (mm  SWE) 
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3  Snow  Depth  and  Density  Modeling 

Snow  depth  and  density  are  simulated  in  GSSHA  independently  of  the 
melting  algorithm  used.  They  are  calculated  hourly  using  routines  taken 
from  the  SNAP  model  (Albert  and  Krajeski  1998)  and  incorporated  within 
GSSHA  as  subroutines.  Information  related  to  SWE,  snow  depth,  density, 
snow  saturation,  effective  porosity,  and  hydraulic  conductivity  are 
exchanged  among  GSSHA  and  the  SNAP  subroutines.  GSSHA  simulates 
the  SWE,  while  the  SNAP  subroutines  simulate  the  remaining  parameters. 
More  information  on  the  connections  between  GSSHA  and  SNAP  and  how 
these  connections  improve  the  vertical  and  lateral  flow  simulation  within 
the  snow  pack  will  be  discussed  in  an  upcoming  report  on  melt-water 
transport  within  GSSHA. 

Albert  and  Krajeski  (1998)  state  the  following  in  regard  to  snow  depth 
calculations  in  the  SNAP  model: 

“The  snow  depth  prediction  equations  are  based  on  the  form  of 
Anderson  (1973),  which  are  also  used  by  Jordan  (1991).  New  snow 
is  calculated  from  precipitation  input  data.  For  snow  settlement  the 
model  predicts  the  rate  of  densification,  calculates  the  overburden 
and  adjusts  the  snow  depth  accordingly.” 

Jordan  (1991)  refers  to  the  SNow  THERMal  (SNTHERM)  model  (Jordan 
1991  and  1998),  often  regarded  as  one  of  the  most  advanced  energy 
balance  snow  models  available. 

The  snow  density  calculated  by  the  SNAP  subroutines  changes  in  response 
to  snow  accumulation,  settlement,  and  melt,  and  is  computed  with  the  snow 
depth  predictions  or  as  updated  by  the  user  (Albert  and  Krajeski  1998). 
GSSHA  currently  does  not  allow  the  user  to  update  the  snow  density  as  a 
calibration  parameter,  leaving  it  to  be  calculated  based  on  the  snow  depth 
predictions  of  the  SNAP  subroutines.  Because  of  this,  the  SNAP  subroutines 
do  not  have  calibrated  parameters  while  running  within  GSSHA. 

For  more  information  on  the  SNAP  model,  consult  Albert  and  Krajeski 
(1998).  For  more  information  on  the  depth  prediction  equations  used 
within  the  SNAP  subroutines,  consult  Anderson  (1973)  and  Jordan  (1991; 
1998). 
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4  Code  Validation 

In  GSSHA,  snow  accumulation  and  melt  are  calculated  on  a  cell-by-cell 
basis.  The  GSSHA  model,  in  relation  to  the  snow  processes,  is  essentially  a 
mosaic  of  point  models.  Therefore,  to  test  and  compare  GSSHA’s  original 
and  new  snow  modeling  capabilities,  several  point  location  models  were 
developed  to  test  the  snow  accumulation  and  melt  algorithms  that  GSSHA 
now  employs  in  each  cell.  By  capturing  snow  accumulation  melt  at  the  point 
scale  accurately,  the  entire  watershed  snow  accumulation  and  melt 
processes  of  the  snowpack  are  tested. 

4.1  Data  collection 

Data  associated  with  nine  sites  in  the  western  United  States  were  collected 
to  support  several  test  cases  wherein  GSSHA  model  simulated  output  from 
the  three  snow  accumulation  and  melt  methods  were  compared  with  their 
observed  counterparts.  Simulated  snow  depth  was  also  compared  with 
observed  snow  depth  data  at  two  sites  at  the  Senator  Beck  Basin  study  area 
in  Colorado. 

Input  meteorological  data  requirements  for  the  EB  and  HY  methods  include 
hourly  precipitation,  temperature,  pressure,  relative  humidity,  and  wind 
speed  data.  The  TI  method  requires  only  hourly  temperature  and 
precipitation  data.  Observed  meteorological  and  snow  data  were  collected 
from  seven  snow  telemetry  (SNOTEL)  sites  (NWCC  2010),  three  study  plots 
within  the  Senator  Beck  Basin  study  area,  and  the  Lake  Tahoe,  California 
airport  (Figure  2).  The  data  collected  from  each  site  were  used  to  compare 
simulated  SWE  results  with  their  observed  counterparts  and  to  force  the 
GSSHA  point  snow  models.  At  two  locations  within  the  Senator  Beck  Basin 
study  area  simulated  snow  depth  was  compared  to  the  observed  data. 

4.1.1  SNOTEL  sites 

SNOTEL  sites  across  the  western  United  States  were  investigated  to 
determine  whether  they  had  the  necessary  data  to  test  the  different  snow 
routines.  Six  SNOTEL  sites  were  selected  for  model  testing.  An  additional 
site  in  the  Trout  Creek  Basin  (TCB),  Hagan  Meadows,  California,  was  also 
chosen  because  it,  along  with  the  HMET  data  from  the  nearby  Lake  Tahoe 
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airport,  had  the  required  data  to  analyze  model-to-measurement  misfit. 
Table  1  shows  the  SNOTEL  site  name,  abbreviation,  elevation,  and 
timeframe  of  the  data  collected  from  each  site. 


Table  1.  SNOTEL  sites. 


Site  Name 

Abbreviation 

Elevation  (m) 

Start  Date 

End  Date 

Sheldon,  NV 

SH 

1786 

10/01/2000 

07/20/2010 

Sierra  Blanca,  NM 

SB 

3133 

10/01/2004 

12/06/2005 

Myrtle  Creek,  ID 

MC 

1073 

10/01/2005 

08/13/2009 

Madison  Butte,  OR 

MB 

1570 

10/01/2005 

07/20/2010 

Lost  Horse,  WA 

LH 

1561 

10/01/2003 

07/20/2010 

Dry  Lake,  CO 

DL 

2560 

10/01/2003 

07/21/2010 

Hagan  Meadows,  CA 

TCB 

2370 

01/01/2000 

09/30/2009 

4.1.2  Lake  Tahoe  airport  data 

Because  the  Hagan’s  Meadow  SNOTEL  site  was  not  equipped  with  the 
correct  sensors  to  measure  the  pressure,  relative  humidity,  wind  speed,  and 
radiation  data,  these  data  were  collected  from  the  Lake  Tahoe  airport 
(Figure  3)  via  the  US  Air  Force  14th  Weather  Squadron.  The  airport  data 
were  collected  at  an  elevation  of  1,911-m  from  January  1,  2000,  to 
September  30,  2009. 
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Figure  3.  Hagan’s  Meadow  and  Lake  Tahoe  airport  locations. 


Although  having  HMET  data  collected  from  two  sites  (KTVL  and  Hagan’s 
Meadow  SNOTEL)  is  not  preferred,  this  is  representative  of  a  typical 
modeling  scenario  in  which  data  are  collected  from  numerous  sources.  The 
distance  between  the  airport  and  the  Hagan’s  Meadow  SNOTEL  site 
(-7.3  km)  and  the  difference  in  elevation  between  the  airport  and  the 
Hagan’s  Meadow  SNOTEL  site  (-0.5  km)  are  both  significant,  but  as 
previously  mentioned,  represent  a  typical  modeling  scenario. 

4.1.3  Senator  Beck  Basin  study  area  data 

The  Senator  Beck  Basin  study  area  (Figure  4)  is  in  the  Ouray  Ranger 
District  of  the  Uncompahgre  National  Forest  in  the  western  San  Juan 
Mountains  of  southwestern  Colorado  (www.snowstudies.org).  Three  study  plots 
are  within  or  near  the  roughly  1  square  mile  basin:  Putney,  Swamp  Angel, 
and  Senator  Beck  Basin  Study  Sites.  Using  hourly  data  from  all  three  sites, 
input  HMET  data  files  were  created  that  represent  the  Swamp  Angel 
(SBSA)  and  Senator  Beck  Basin  Study  Site  (SBSB).  Staff  at  the  Center  for 
Snow  and  Avalanche  Studies  approximately  15  times  per  year  measured 
the  SWE,  snow  depth  (SD),  and  snow  density  at  SBSA  and  SBSB.  Table  2 
shows  which  study  plot  datasets  were  used  to  create  the  SBSB  and  SBSA 
datasets.  Senator  Beck  Basin  study  area  data  was  collected  from  October,  1 
2006,  through  June  24,  2010. 
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Table  2.  SBSB  and  SBSA  datasets  and  data  sources. 


SBSB  Dataset 

SBSA  Dataset 

Temperature 

B 

C 

Precipitation 

c 

C 

Relative  Humidity 

A 

A 

Pressure 

c 

c 

Wind  Speed 

A 

A 

Snow  Water  Equivalent 

B 

C 

Snow  Depth 

B 

c 

Snow  Density 

B 

C 

A  -  Data  came  from  Putney  Study  Plot  (3756-m) 

B  -  Data  came  from  Senator  Beck  Basin  Study  Plot  (3714-m) 
C  -  Data  came  from  Swamp  Angel  Study  Plot  (3371-m) 


4.2  Metrics 

Accuracy  among  the  snow  melt  models  is  quantified  by  comparing  the 
Nash-Sutcliffe  efficiency  (NSE,  Equation  23)  (Nash  and  Sutcliffe  1970), 
root  mean  square  error  (RMSE,  Equation  24),  and  mean  absolute  error 
(MAE,  Equation  25).  The  observed  and  simulated  data  were  compared  for 
the  entirety  of  the  simulation,  and  for  all  snow  events  combined 
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sequentially.  One  snow  event  was  defined  per  water  year  (01  Oct  through 
31  Sept),  and  was  defined  as  the  largest  period  with  continuously  observed 
snow.  Analysis  of  snow  events  was  carried  out  to  not  bias  the  results  with 
numerous  snowless  days. 


JVSE  =  1-^1^°  ^4 

EJqs-Qo) 


RMSE  =  - 


T- 1 


mae=±±Iq^-q: 

1  t= 1 


(23) 


(24) 


(25) 


where 

NSE  =  efficiency  of  simulation 
RMSE  =  root  mean  square  error  of  simulation 
MAE  =  mean  absolute  error  of  simulation 
T  =  time  span  of  period  of  interest 
t  =  time  step  within  time  span,  T 
Qq  =  observed  value  at  time  t 
Qm  =  modeled  value  at  time  t 
Q0  =  average  observed  value  over  time  span,  T 

4.3  Calibration  of  the  snow  models 

The  EB  model  requires  two  total  calibration  parameters  -  the  snow  adjust¬ 
ment  factor  ( SCF ),  which  accounts  for  misreadings  in  precipitation  of  the 
gauge  system,  and  the  temperature  at  which  precipitation  begins  to  fall  in 
the  form  of  snow  (. MBASE ).  The  HY  model  requires  four  total  calibrated 
parameters,  two  associated  with  the  snow  accumulation  ( SCF  and  MBASE), 
and  two  associated  with  the  heat  deficit  ( TIPM  and  NMf).  The  TI  model 
requires  nine  total  calibrated  parameters  associated  with  snow  accumula¬ 
tion,  heat  deficit,  melt,  and  water  distribution  within  the  snowpack.  These 
calibrated  parameters  are  in  Section  2.2.  All  three  snow  models  were 
calibrated  using  observed  SWE  data.  The  PEST  model  (Doherty  et  al.  2002) 
was  used  to  calibrate  the  snow  models,  and  the  lower  and  upper  limits  used 
for  each  parameter  are  in  Table  3.  The  calibrated  parameters  for  the  EB, 
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HY,  and  TI  snow  models  are  in  Tables  4,  5,  and  6,  respectively.  Based  on  the 
PEST  calibration  of  forward  model  calls,  the  EB  method  is  most  sensitive  to 
SCF ;  the  HY  method  is  most  sensitive  to  TIPM  and  NMf,  and  the  TI  method 
is  most  sensitive  to  MFmin  and  SCF. 


Table  3.  Calibration  parameter  limits  and  units. 


Calibration 

Parameter 

Description 

Lower  Limit 

Upper  Limit 

Units 

MFmax 

maximum  melt  factor 

0.00 1 

2.50  1 

(mm^C1-  6-hr1) 

MFmin 

minimum  melt  factor 

0.00 1 

1.00  1 

(mm^C1-  6-hr1) 

ffuse 

fraction  of  precipitation  as  rain 

0.00° 

1.00° 

unitless 

SCF 

snow  adjustment  factor 

0.80° 

1.20° 

unitless 

TIPM 

weighting  factor  for  heat  deficit 

0.01 1 

1.00  1 

unitless 

NMf 

negative  melt  factor 

0.00° 

1.00° 

(mm^C-MS-hr1). 

fua 

average  wind  function 

0.00° 

1.00° 

(mm-mb^e-hr1) 

PLWHC 

percent  liquid  water  holding 
capacity 

0.00  1 

0.40  1 

unitless 

MBASE 

temperature  at  which 
precipitation  begins  to  fall  as 
snow 

0.00  2 

5.00  2 

°c 

°Physical  limits  of  parameters  or  based  on  previous  calibration  exercises 
^-Snow  accumulation  and  ablation  model  -  Snow- 17  (Anderson  2006) 
2National  snow  and  ice  data  center  website  (NSIDC  2012) 


Table  4.  EB  snow  model  calibration  parameters. 


DL 

LH 

MB 

MC 

SB 

SH 

TCB 

SBSA 

SBSB 

SCF 

1.20 

1.20 

1.20 

1.20 

1.20 

1.20 

1.20 

1.20 

1.04 

MBASE 

4.02 

4.05 

4.11 

4.82 

3.79 

4.12 

4.00 

4.00 

4.00 

*most  sensitive  parameter  highlighted 

DL  -  Dry  Lake,  CO  LH  -  Lost  Horse,  WA  MB  -  Madison  Butte,  OR  MC  -  Myrtle  Creek,  ID  SB  -  Sierra  Blanca,  NM  SH  -  Sheldon,  NV 
TCB  -  Hagan's  Meadow,  CA  SBSA  -  Swamp  Angel  Study  Site,  CO  SBSB  -  Senator  Beck  Basin  Study  Site,  CO 


Table  5.  HY  snow  model  calibration  parameters. 


DL 

LH 

MB 

MC 

SB 

SH 

TCB 

SBSA 

SBSB 

SCF 

1.07 

1.19 

1.20 

0.80 

1.03 

0.88 

0.98 

1.20 

0.83 

TIPM 

0.00 

0.00 

0.21 

0.01 

0.29 

0.10 

0.00 

0.30 

0.02 

NMf 

0.17 

0.17 

0.55 

0.53 

0.42 

0.37 

0.58 

0.28 

0.18 

MBASE 

3.98 

4.12 

3.17 

4.23 

3.74 

3.80 

4.00 

4.00 

4.00 

*two  most  sensitive  parameters  highlighted 
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Table  6.  Tl  snow  model  calibration  parameters. 


DL 

LH 

MB 

MC 

SB 

SH 

TCB 

SBSA 

SBSB 

MFmax 

0.56 

0.60 

0.57 

0.52 

0.57 

0.12 

0.70 

0.55 

0.77 

MFmin 

0.01 

0.00 

0.01 

0.02 

0.00 

0.10 

0.00 

0.01 

0.00 

fruse 

0.00 

0.00 

0.02 

0.00 

0.03 

0.28 

0.01 

0.01 

0.00 

SCF 

1.00 

1.00 

1.02 

0.94 

1.01 

0.80 

0.92 

1.00 

0.80 

TIPM 

0.00 

0.36 

0.00 

0.00 

0.28 

1.00 

0.00 

0.20 

0.24 

NMf 

0.18 

0.34 

0.25 

0.05 

0.23 

0.00 

1.00 

0.20 

0.26 

fua 

0.90 

0.90 

0.90 

0.90 

0.90 

0.11 

0.90 

0.18 

0.11 

PLWHC 

0.16 

0.29 

0.22 

0.00 

0.15 

0.01 

0.40 

0.20 

0.10 

MBASE 

4.00 

3.97 

3.98 

4.05 

3.95 

4.27 

4.00 

4.00 

4.00 

*three  most  sensitive  parameters  highlighted 
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5  Results 

The  EB,  HY,  and  TI  snow  models  were  evaluated  using  observed  SWE  data. 
As  discussed  in  Section  4.2,  the  simulations  were  compared  for  the  entire 
simulation  (Overall)  and  each  major  snow  event  for  each  water  year 
(Events).  Because  SBSA  and  SBSB  datasets  have  periodic  observations 
(approximately  15  measurements  per  year)  instead  of  hourly  observations, 
each  was  considered  an  event,  as  reflected  in  the  Events  statistics.  Because 
of  the  few  time  steps  with  observation  data,  the  Overall  statistics  were  not 
calculated  for  the  SBSA  and  SBSB  datasets.  The  modeling  fit  statistics  from 
all  the  tests  are  listed  in  Tables  7  through  15.  The  closer  the  efficiency  value 
is  to  1.0  and  the  closer  RMSE  (m)  and  MAE  (m)  are  to  0.0,  the  more 
accurate  the  simulation.  Figures  5  through  12  and  Figure  14  show  the 
observed  and  simulated  SWE  from  all  three  snow  models.  Figures  13  and  15 
show  the  observed  and  simulated  SD  from  all  three  snow  models  for  the 
SBSA  and  SBSB  sites. 


Table  7.  Dry  Lake,  Colorado  site  SWE  fit  statistics. 


Overall 

Event 

Model 

NSE 

RMSE 

MAE 

NSE 

RMSE 

MAE 

EB 

0.83 

0.09 

0.05 

0.72 

0.11 

0.08 

HY 

0.93 

0.06 

0.03 

0.89 

0.07 

0.05 

TI 

0.95 

0.05 

0.03 

0.92 

0.06 

0.04 

^Shaded  indicate  the  best  fit. 
NSE  (%),  RMSE  (m),  and  MAE  (m) 


Table  8.  Lost  Horse,  Washington  site  SWE  fit  statistics. 


Overall 

Event 

Model 

NSE 

RMSE 

MAE 

NSE 

RMSE 

MAE 

EB 

0.76 

0.08 

0.05 

0.60 

0.11 

0.08 

HY 

0.86 

0.06 

0.04 

0.82 

0.08 

0.05 

TI 

0.89 

0.06 

0.03 

0.87 

0.06 

0.04 

^Shaded  indicate  the  best  fit. 
NSE  (%),  RMSE  (m),  and  MAE  (m) 


ERDC/CHL  TR-13-4 


21 


Table  9.  Madison  Butte,  Oregon  site  SWE  fit  statistics. 


Overall 

Event 

Model 

NSE 

RMSE 

MAE 

NSE 

RMSE 

MAE 

EB 

0.05 

0.06 

0.03 

-1.28 

0.11 

0.08 

HY 

0.21 

0.06 

0.03 

-0.47 

0.09 

0.06 

Tl 

0.44 

0.05 

0.02 

0.63 

0.04 

0.04 

^Shaded  indicate  the  best  fit. 
NSE  (%),  RMSE  (m),  and  MAE  (m) 


Table  10.  Myrtle  Creek,  Idaho  site  SWE  fit  statistics. 


Overall 

Event 

Model 

NSE 

RMSE 

MAE 

NSE 

RMSE 

MAE 

EB 

0.39 

0.03 

0.01 

-0.05 

0.05 

0.03 

HY 

0.93 

0.01 

0.00 

0.89 

0.02 

0.01 

Tl 

-4.13 

0.07 

0.04 

-4.44 

0.08 

0.05 

^Shaded  indicate  the  best  fit. 
NSE  (%),  RMSE  (m),  and  MAE  (m) 


Table  11.  Sierra  Blanca,  New  Mexico  site  SWE  fit  statistics. 


Overall 

Event 

Model 

NSE 

RMSE 

MAE 

NSE 

RMSE 

MAE 

EB 

0.66 

0.07 

0.03 

0.35 

0.11 

0.08 

HY 

0.75 

0.06 

0.02 

0.54 

0.10 

0.07 

Tl 

0.79 

0.05 

0.02 

0.65 

0.08 

0.05 

^Shaded  indicate  the  best  fit. 
NSE  (%),  RMSE  (m),  and  MAE  (m) 


Table  12.  Sheldon,  Nevada  site  SWE  fit  statistics. 


Overall 

Event 

Model 

NSE 

RMSE 

MAE 

NSE 

RMSE 

MAE 

EB 

0.25 

0.01 

0.00 

-1.93 

0.02 

0.01 

HY 

-0.76 

0.01 

0.00 

-0.25 

0.01 

0.01 

Tl 

0.18 

0.01 

0.00 

0.19 

0.01 

0.01 

^Shaded  indicate  the  best  fit. 
NSE  (%),  RMSE  (m),  and  MAE  (m) 
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Table  13.  Hagan’s  Meadow,  California  site  SWE  fit  statistics. 


Overall 

Event 

Model 

NSE 

RMSE 

MAE 

NSE 

RMSE 

MAE 

EB 

0.37 

0.13 

0.06 

-0.35 

0.19 

0.14 

HY 

0.78 

0.07 

0.03 

0.75 

0.08 

0.05 

TI 

0.82 

0.07 

0.04 

0.67 

0.09 

0.07 

^Shaded  indicate  the  best  fit. 
NSE  (%),  RMSE  (m),  and  MAE  (m) 


Table  14.  Swamp  Angel,  Colorado  site  SWE  fit  statistics. 


Overall 

Event 

Model 

NSE 

RMSE 

MAE 

NSE 

RMSE 

MAE 

EB 

N/A 

N/A 

N/A 

-1.95 

0.39 

0.31 

HY 

N/A 

N/A 

N/A 

0.10 

0.21 

0.11 

TI 

N/A 

N/A 

N/A 

0.82 

0.10 

0.07 

^Shaded  indicate  the  best  fit. 
NSE  (%),  RMSE  (m),  and  MAE  (m) 


Table  15.  Senator  Beck  Basin,  Colorado  site  SWE  fit  statistics. 


Overall 

Event 

Model 

NSE 

RMSE 

MAE 

NSE 

RMSE 

MAE 

EB 

N/A 

N/A 

N/A 

0.41 

0.16 

0.14 

HY 

N/A 

N/A 

N/A 

0.70 

0.11 

0.08 

TI 

N/A 

N/A 

N/A 

0.61 

0.13 

0.11 

^Shaded  indicate  the  best  fit. 
NSE  (%),  RMSE  (m),  and  MAE  (m) 


5.1  Dry  Lake,  Colorado  site  results 

As  shown  in  Table  7,  this  site  had  the  highest  fit  between  observed  and 
simulated  SWE.  It  also  had  one  of  the  highest  SWE  volumes  with  water  year 
SWE  maximums  ranging  from  0.4  m  and  almost  1.0  m.  As  expected,  Figure 
5  shows  that  the  EB  method  uniformly  underestimated  SWE,  while  the  HY 
and  TI  methods  typically  followed  well  with  the  observed  SWE.  The 
calibrated  snow  adjustment  factor  ( SCF )  for  the  HY  and  TI  methods  was 
around  1.0,  while  the  EB  method  had  a  value  of  1.2.  This  highlights  the  fact 
that  the  EB  method  overestimates  melt  even  when  a  high  SCF  value  is  used. 
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Figure  5.  Dry  Lake,  Colorado  SWE  results  (B). 
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Figure  6.  Lost  Horse,  Washington  SWE  results  (B). 


5.2  Lost  Horse,  Washington  site  results 

The  HY  and  TI  simulated  SWE  matched  well  with  the  observed  SWE  at  the 
Lost  Horse  site.  Even  with  a  good  fit,  the  HY  and  TI  methods  did  prolong 
the  snow  past  observed  during  Events  l  and  4  through  9  (Figure  6).  It  is 
evident,  at  this  site  in  particular,  how  the  EB  underestimates  the  SWE 
compared  to  the  HY  and  TI  methods,  as  evidenced  in  Figure  6  and  Table  8 
(EB  has  the  lowest  NSE  values  and  the  highest  RMSE  and  MAE  values). 
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Figure  7.  Madison  Butte,  Oregon  SWE  results. 


Figure  8.  Myrtle  Creek,  Idaho  SWE  results. 
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Figure  9.  Sierra  Blanca,  New  Mexico  SWE  results. 


5.3  Madison  Butte,  Oregon,  site  results 

The  overall  results  from  this  site  were  poor.  All  three  methods 
underestimated  the  SWE  and  the  snowpack  was  simulated  to  last  longer 
than  observed  (Figure  7).  The  relatively  small  snowpack  (SWE  maximum  of 
0.3  m)  might  have  skewed  the  results  (even  a  relatively  small  error  is 
magnified  in  such  a  small  snowpack).  As  evidenced  in  Table  9,  the  TI 
method  outperformed  the  energy  balance  based  methods.  This  might  be  due 
to  the  ability  of  the  TI  method  to  better  tune  the  model  to  observed  data 
because  it  has  nine  calibrated  parameters,  rather  than  the  four  of  the  HY 
method  and  two  of  the  EB  method. 

5.4  Myrtle  Creek,  Idaho  site  results 

The  calibrated  TI  method  overestimated,  the  EB  underestimated,  and  the 
HY  method  matched  the  SWE  very  well  during  every  event  at  this  site.  The 
snow  volumes  were  small  (SWE  max  of  0.15  m),  but  the  HY  method 
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performed  well,  as  seen  in  Figure  8  and  Table  10.  The  poor  results  of  the  TI 
method  could  indicate  an  issue  in  basing  all  melt  processes  on  one 
parameter,  temperature.  More  investigation  into  this  site  might  be  required 
to  fully  determine  the  cause  of  this  poor  performance  by  the  TI  method. 


Figure  10.  Sheldon,  Nevada  SWE  results  (A). 
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Figure  11.  Hagan’s  Meadows,  California  SWE  results  (A). 
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Figure  12.  Swamp  Angel  study  plot,  Colorado,  SWE  results. 


Figure  13.  Swamp  Angel  study  plot,  Colorado,  snow  depth  results. 
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Figure  14.  Senator  Beck  Basin  study  plot,  Colorado,  SWE  results. 


Figure  15.  Senator  Beck  Basin  study  plot,  Colorado,  snow  depth  results. 
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5.5  Sierra  Blanca,  New  Mexico,  site  results 

The  first  event  (Figure  9)  at  this  site  was  best  simulated  using  the  HY 
method.  All  other  events  were  modeled  best  using  the  TI  method. 
Compared  to  the  other  sites,  the  Sierra  Blanca  site  had  a  medium 
snowpack  with  a  maximum  SWE  of  0.5  m. 

5.6  Sheldon,  Nevada,  site  results 

This  site  had  the  worst  results  of  all  sites  analyzed.  It  also  had  the  smallest 
snowpack  of  any  site,  with  a  maximum  SWE  of  just  0.042-m.  The  EB 
method  performed  best  in  Events  4  through  6,  (Figure  10).  This  site 
indicates  the  EB  method  might  perform  well  in  shallow  snowpacks  where 
overcoming  the  heat  deficit  prior  to  melt  is  not  physically  relevant. 

5.7  Hagan’s  Meadow,  California,  site  results 

This  site  had  the  longest  period  of  record  of  any  of  the  test  sites.  As 
evidenced  in  Figure  11  and  Table  13,  both  the  HY  and  TI  methods 
performed  consistently  well  at  this  site.  The  EB  method  consistently 
underestimated  the  SWE.  The  HY  method  accurately  captured  the 
accumulation  and  melt  in  Events  4,  7,  and  8,  and  over-  or  under¬ 
estimated  the  SWE  in  the  other  events.  In  Events  3,  5,  and  9  the  HY 
method  prolonged  the  snowpack  much  longer  than  the  observed  record. 
This  prolonging  of  the  snowpack  is  not  captured  within  the  Event  statistics 
because  they  are  only  calculated  while  snow  is  present  in  the  observed 
record.  For  this  reason,  despite  the  Event  results  in  Table  13,  the  TI 
method  likely  outperformed  the  HY  method. 

5.8  Swamp  Angel,  Colorado,  site  results 

For  the  Swamp  Angel  site,  the  TI  method  performed  the  best,  followed  by 
the  HY  and  the  EB  methods.  Sporadic  observed  SWE  readings  prevented 
the  Overall  statistics  from  being  computed.  The  Event  fit  statistics  (Table 
14)  show  that  the  TI  method  performed  best.  Visually  from  Figure  12,  the 
HY  method  performs  very  well  except  for  the  melting  in  the  2007  to  2008 
plot.  Although  the  SWE  results  were  accurate,  Figure  13  shows  the  snow 
depth  was  consistently  underestimated. 
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5.9  Senator  Beck  Basin,  Colorado,  site  results 

For  this  site,  the  HY  method  performed  the  best,  followed  by  the  77  and 
EB  methods.  Sporadic  observed  SWE  readings  prevented  the  Overall 
statistics  from  being  computed.  The  Event  fit  statistics  (Table  15)  show  the 
HY  method  performs  best.  Visually  from  Figure  14,  all  methods  slightly 
overestimated  SWE,  except  for  the  2009  to  2010  plot  when  SWE  was 
underestimated.  Figure  15  shows  that  the  snow  depth  results  were 
accurate,  except  for  the  2009  to  2010  plot  where  the  simulations 
underestimated  snow  depth. 
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6  Conclusions 

The  results  presented  clearly  demonstrate  the  new  snow  models, 
temperature-index  (TI)  and  hybrid  energy  balance  (HY),  now  present 
within  GSSHA  outperform  the  original  snow  model,  energy  balance  (EB), 
which  consistently  underestimated  SWE.  With  the  addition  of  the  heat 
deficit,  the  snowpack  tends  to  accumulate  more  and  melt  slower  than  when 
the  original  EB  method  was  used,  in  most  cases  making  the  simulation 
more  accurate.  The  importance  of  including  the  heat  deficit  is  further 
evidenced  by  the  EB  model  performing  best  at  sites  with  a  small  snowpack. 
At  these  sites  accounting  for  the  heat  deficit  will  make  a  minimal  difference 
because  the  small  snowpack  likely  has  a  small  heat  deficit  to  overcome.  In 
comparison,  sites  with  larger  snowpacks,  and  therefore  a  larger  heat  deficit, 
were  best  simulated  by  the  TI  and  HY  models,  which  account  for  the  heat 
deficit. 

For  most  of  the  simulations,  the  TI  model  performed  best,  but  this  model 
requires  nine  calibration  parameters.  When  deploying  a  model  in  a  data- 
poor  environment,  use  of  the  HY  model  is  advised  because  it  performs  well 
with  only  four  calibration  parameters  (included  in  the  four  are  two 
associated  with  snow  accumulation  routines,  which  are  not  specific  to  the 
melt  method  deployed).  At  sites  with  a  small  snowpack  (Sheldon,  Nevada 
and  Myrtle  Creek,  Idaho)  the  EB  and  HY  models  outperformed  or 
performed  nearly  as  well  as  the  TI  model.  For  this  reason,  use  of  the  HY 
model  is  advised  in  areas  with  a  small  maximum  snowpack  (-0.15-01  SWE 
or  less). 

The  snow  depth  simulations  shown  in  Figures  13  and  15  show  that  when 
the  SWE  is  simulated  accurately,  the  SNAP  model  (in  the  form  of 
subroutines  now  within  GSSHA)  simulates  the  snow  depth  reasonably. 
The  results  show  that  the  SNAP  subroutines  might  slightly  overestimate 
snow  depth  when  compared  to  observed  data. 

Overall,  the  GSSHA  model  has  been  improved  by  the  inclusion  of  the  HY 
and  TI  snow  melt  methods,  as  well  as  modifications  to  how  snow  fall 
accumulation  is  simulated.  GSSHA  can  now  more  accurately  simulate 
snow  accumulation  and  melt  within  each  cell,  better  capturing  the 
snowpack  and  melt  over  an  entire  watershed. 
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